中亚地区陆表物候逐年数据集（1982–2015）
马勇刚，刘素红*
新疆大学资源与环境科学学院，乌鲁木齐
830046
摘 要：中亚（30°N–60°N，50°E–100°E）为世界典型的干旱半干旱地区，研究该区域的植物物候变化对理解水分胁迫环境下地表植被对气候变化的响应具有重要意义。本数据集在GIMMS NDVI3g.v1长时间序列NDVI数据的基础上，利用阈值法和拐点法两种遥感物候提取方法对中亚地区1982–2015年共34年的地表植物物候信息进行提取。该数据集空间参考系统为地理坐标系，空间分辨率为0.083,3°，约8 km。数据集含有两个子集：在预先对时间序列进行逻辑函数拟合进行数据重构后，V1数据集为采用20%动态阈值法计算得到生长季的开始期、停止期和生长期长度，V2数据集为在采用拐点法计算得到生长季的开始期、停止期和生长期长度。V1数据集包含204个文件，数据格式为.hdr（头文件）和.img（数据文件）；V2数据集包含204个文件，数据格式为.hdr（头文件）和.tif（数据）。本数据集包含2个文件夹共408个文件，数据量为126 MB（压缩为一个文件，数据量为75.6 MB）。
关键词：陆表物候；遥感；中亚
DOI: 10.3974/geodp.2020.01.05
1 前言
全球变暖对地球生态系统的影响是本世纪地球系统科学领域被广泛关注的议题之一。植物物候信息作为佐证气候变化的重要参考被广泛的收集、归纳和分析，以此判断历史气候变化的特征和趋势[1]。与此同时，地球系统模拟者也急需要准确的植物物候空间数据和其内在机理参与地球系统模拟运算，从而更精确的模拟区域尺度下物候在当前乃至未来长期的气候变化过程中所表现出的响应，并以此为中间参量实现大区域甚至是地球系统能量物质运移的科学模拟和计算[2–3]。
相对于巨大需求而言，陆表物候的数据来源相对缺乏。目前植物物候数据来源可以分为人工定点观测、遥感数据反演以及基于数字相机或者无人机原位观测的物候反演，三种方法各有优劣。人工定点观测能够准确获取植物物种尺度的物候信息，但观测范围很小，很难形成大区域尺度的同步观测，且易受人为因素的影响；遥感数据的反演是利用可见光至近红外植物生长阶段的光谱响应特征，构建时间尺度上的植物生长曲线，通过提取曲线特征值确定植物物候数据，这种方法能够形成景观尺度或区域尺度的物候信息数据，但与人为观测数据存在尺度差异，需要进行仔细验证[4]；数字相机和无人机等利用现代技术手段进行物候测量，其原理与遥感物候反演类似，虽然时间和空间尺度的分辨率较高，但缺乏长时间序列的支持，主要还用于对替代人为观测物候数据、与遥感物候进行交叉验证及多尺度分析工作[5–7]。
中亚地区是世界最主要的干旱半干旱地区之一，在这一地区开展植物物候研究工作能够加深我们理解水胁迫条件下植物物候对气候变化的响应机制。目前在中亚地区除极少数农业气象站和植物园开展面向特定植物种物候观测工作外，其他人为观测数据很难获取。遥感物候数据仍是中亚地区物候数据的主要来源。GIMMS NDVI3g.v1数据是目前全球发布的时间跨度最长的植被指数数据，利用该数据进行物候数据的反演能够从30年尺度上进行本区域范围内物候数据的科学分析。因此，本数据集采用目前被广泛应用的阈值法和拐点法反演了中亚地区1982–2015年的开始期、停止期和生长季长度等陆表物候数据集。
2 数据集元数据简介
《中亚地区陆表物候逐年数据集（1982–2015）》[8]的主要元数据要素见表1。
表1
《中亚地区陆表物候逐年数据集（1982–2015）》元数据简表
|
条 目 |
描 述 |
|
数据集名称 |
中亚地区陆表物候逐年数据集（1982–2015） |
|
数据集短名 |
LSP_CA |
|
作者信息 |
马勇刚 AAH-5322-2019, 新疆大学资源与环境科学学院, mayg@xju.edu.cn |
|
|
刘素红 AAB-2538-2020, 新疆大学资源与环境科学学院, liush@bnu.edu.cn |
|
地理区域 |
30°N–60°N，50°E–100°E |
|
数据年代 |
1982–2015 |
|
时间分辨率 |
年 |
|
空间分辨率 |
0.083,3°（约8 km） |
|
数据格式 |
.img、.tif、.hdr |
|
数据量 |
126 MB（压缩后75.6 MB） |
|
数据集组成 |
2个文件夹下共计408个数据文件压缩为1个文件 （1）V1文件夹下：阈值法获取开始期、停止期和生长季长度数据，共204个文件，其中102个.img实体数据，102个.hdr头文件 （2）V2文件夹下：拐点法获取开始期、停止期和生长季长度数据，共204个文件，其中102个.tif实体数据，102个.hdr头文件 |
|
基金项目 |
中华人民共和国科学技术部（2017YEF0118100）；国家自然科学基金（41761013，41861053）；新疆自治区高校计划项目（XJEDU2017M007） |
|
数据计算环境 |
Matlab 2014b |
|
出版与共享服务平台 |
全球变化科学研究数据出版系统 http://www.geodoi.ac.cn |
|
地址 |
北京市朝阳区大屯路甲11号100101，中国科学院地理科学与资源研究所 |
|
数据共享政策 |
全球变化科学研究数据出版系统的“数据”包括元数据（中英文）、实体数据（中英文）和通过《全球变化数据学报》（中英文）发表的数据论文。其共享政策如下：（1）“数据”以最便利的方式通过互联网系统免费向全社会开放，用户免费浏览、免费下载；（2）最终用户使用“数据”需要按照引用格式在参考文献或适当的位置标注数据来源；（3）增值服务用户或以任何形式散发和传播（包括通过计算机服务器）“数据”的用户需要与《全球变化数据学报》（中英文）编辑部签署书面协议，获得许可；（4）摘取“数据”中的部分记录创作新数据的作者需要遵循10%引用原则，即从本数据集中摘取的数据记录少于新数据集总记录量的10%，同时需要对摘取的数据记录标注数据来源[9] |
|
数据和论文检索系统 |
DOI，DCI，CSCD，WDS/ISC，GEOSS，China GEOSS，Crossref |
3 数据研发方法
3.1 算法原理
对GIMMS
NDVI3g.v1数据首先进行逻辑函数拟合和数据重构。
(1)
式中，NDVImax为NDVI最大值，NDVImin为NDVI最小值，t为儒略日(d)，VI(t)是拟合后t日的NDVI值，A和B为拟合参数。NDVImax与NDVImin之差为植被生长曲线振幅。
（1）动态阈值法[10]
对拟合后数据进行判断，当儒略日t的NDVI值大于等于20%的振幅时，则认定t为开始期时间，当儒略日t的NDVI值小于等于20%的振幅时，则认定t为停止期时间。
(2)
（2）拐点法[11]
对拟合后数据时序计算曲率变化率。曲率计算如下:
(3)
其中
(4)
曲率变化率为:
(5)
通过计算曲率变化率
的局部极值，即可确定拐点位置，从而确定开始期和停止期。
3.2 技术路线
数据集研发的主要过程是：对GIMMS NDVI3g.v1数据进行相应的数据预处理，获取中亚地区1982–2015年NDVI时间序列数据，继而通过逻辑函数拟合获取各栅格的植物生长拟合曲线，利用阈值法和拐点法进行了2种物候数据集的提取，并得到2套包含开始期、停止期和生长季长度3种物候参数的子数据集，最终通过数据集成，完成LSP_CA中亚地区的陆表物候数据集的制作。研发技术路线示于图1。
4 数据结果与验证
4.1
数据集组成
《中亚地区陆表物候逐年数据集（1982–2015）》由阈值法和拐点法2种物候提取方法获取的开始期、停止期和生长季长度数据组成，数据集组成文件的命名方式、数据描述、数据格式、文件个数及数据量见表2。
图1 中亚地区陆表物候逐年数据集（1982–2015）研发技术路线图
表2
中亚地区陆表物候逐年数据集（1982–2015）组成文件简表
|
组成文件 |
命名方式 |
数据描述 |
数据格式 |
文件数 |
数据量 |
|
头文件 |
与物候栅格数据一致 |
行列数，数据类型，空间参考系统 |
.hdr |
204 |
132.3 KB |
|
物候栅格数据 |
LOG_20%_SOSyyyy.img |
阈值法获取的开始期 |
.img |
34 |
28 MB |
|
LOG_20%_EOSyyyy.img |
阈值法获取的停止期 |
.img |
34 |
28 MB |
|
|
LOG_20%_LOSyyyy.img |
阈值法获取的生长期长度 |
.img |
34 |
28 MB |
|
|
LOG_inflexion_SOSyyyy.tif |
拐点法获取的开始期 |
.tif |
34 |
14.1 MB |
|
|
LOG_inflexion_EOSyyyy.tif |
拐点法获取的停止期 |
.tif |
34 |
14.1 MB |
|
|
LOG_inflexion_LOSyyyy.tif |
拐点法获取的生长期长度 |
.tif |
34 |
14.1 MB |
（1）数据头文件（.hdr）。内含所有对应栅格数据的头文件信息，包括数据类型，行列数，投影信息。
（2）物候栅格数据。包括开始期、停止期和生长季长度3种物候参数。以阈值法获取的开始期、停止期和生长期长度栅格数据分别为LOG_20%_SOSyyyy.img、LOG_20%_ EOSyyyy.img和LOG_20%_LOSyyyy.img；以拐点法获取的开始期、停止期和生长期长度栅格数据分别为LOG_inflexion_SOSyyyy.tif、LOG_inflexion_EOSyyyy.tif和LOG_inflexion_ LOSyyyy.tif。其中，yyyy代表四位数年，每一物候栅格数据均有对应名称的数据头文件（.hdr），可在ENVI软件进行导入导出，栅格数据值代表当年的对应物候参数儒略日时间，无效值为0，空间坐标系为经纬度地理坐标。
4.2 数据结果
阈值法和拐点法获得的物候参数的部分示意图示于图2。可见，由阈值法提取的开始
图2 中亚地区陆表物候逐年（1982–2015）数据可视化图
期，停止期和生长季长度等参数在空间上反映了整个中亚区域的空间分布情况，但阈值法相较于拐点法的结果更加平滑。这主要是由NDVI数据在生长曲线上升初期波动较为明显，继而导致曲率变化率异质性较大，空间平滑度较差。
两种方法获得的数据结果显示，中亚范围内34年来的开始期、停止期、生长季长度约73%–85%的范围变化均不显著，但均检测到有15%的区域开始期有显著提前的趋势；两种方法对停止期的结果分析不一致，继而造成对生长季长度的变化分析截然相反。阈值法结果表明，7.3%的区域生长季长度缩短，18.8%的区域生长季长度延长；拐点法则表明，超过21%的区域生长季长度延长，只有3.5%的区域生长季长度缩短。
目前中亚地区可获取的物候数据较少，公开发布的数据为MODIS数据研发的全球MCD12Q2的005和006两个版本的数据集，005版本主要采用拐点法进行数据计算，006版本则主要采用阈值法计算。与MCD12Q2数据相比，本数据集拥有更长的时间分辨率，更适用于大范围均值地表区域的长时间变化分析；但由于本数据空间分辨率较粗，较之MCD12Q2数据更易受到地表异质性的影响，在地表覆盖类型多样或人类干扰强度大的区域可能产生较大误差。
4.3 数据结果分析
物候数据的时间变化是物候研究的重点。为了分析本区域两种数据结果的时间趋势特征，采用Mann-kendall趋势检验[12]，对本数据集34年来两种方法获取的物候数据进行检验，统计结果如表3。
表3 中亚地区陆表物候34年（1982–2015）趋势变化表（%）
|
趋势类型 |
阈值法 |
拐点法 |
||||
|
开始期 |
停止期 |
生长季长度 |
开始期 |
停止期 |
生长季长度 |
|
|
极显著提前（缩短） |
9.65 |
4.47 |
3.43 |
6.81 |
1.77 |
11.38 |
|
显著提前（缩短） |
7.80 |
4.43 |
3.91 |
9.22 |
1.39 |
9.79 |
|
不显著 |
74.50 |
77.95 |
73.86 |
80.47 |
85.68 |
75.36 |
|
显著推迟（延长） |
3.80 |
6.45 |
8.12 |
2.19 |
5.93 |
2.23 |
|
极显著提前（延长） |
4.25 |
6.71 |
10.68 |
1.32 |
3.84 |
1.24 |
5 讨论和总结
本数据集研发的34年中亚地区陆表物候数据库有利于辅助开展本地区气候变化和陆表植物生态系统监测的研究工作。基于GIMMS NDVI3g.v1的物候参数集由于空间分辨率（8 km）较大，必然受到空间异质性问题的干扰，该数据的精度暂时还无法从人工观测的角度进行准确评估。通过实地考察，在中亚地区选择一些具有较高地表一致性且地面覆盖类型长期未发生变化的区域进行实地观测数据的采集，将会对本数据集的精度提供更加准确的评估。此外，使用其他物候数据集如MCD12Q2产品，或基于一些长时间中分辨率遥感时间序列数据如Landsat及哨兵等开展新的物候产品数据集研发，深入尺度转换和分析研究是在该领域下一步的研究工作。本数据集能够从空间变化的角度上揭示一些物候在时空尺度上的变化特征，但容易受到地表覆盖类型变化的影响和干扰。在使用本数据做进一步分析之前，利用历史土地利用覆被数据如MCD12Q1等数据进行影像剔除，将有利于提升研究结论的可靠性。
作者分工：马勇刚对数据集的开发做了总体设计，采集和处理了GIMMS NDVI3g.v1数据；刘素红设计了模型和算法，并做了数据验证；马勇刚撰写了数据论文；刘素红审核了数据论文。
参考文献
[1]
Chuine, I.,
Yiou, P., Viovy, N., et al.
Historical phenology: grape ripening as a past climate indicator [J]. Nature, 2004, 432(7015): 289–290.
[2]
Migliavacca,
M., Sonnentag, O., Keenan, T. F., et al. On the uncertainty of phenological responses
to climate change, and implications for a terrestrial biosphere model [J]. Biogeosciences, 2012, 9(6): 2063.
[3]
Randerson,
J. T., Hoffman, F. M., Thornton, P. E., et
al. Systematic assessment of terrestrial biogeochemistry in coupled climate-carbon
models [J]. Global Change Biology,
2009, 15(10): 2462–2484.
[4]
Cleland, E.
E., Chuine, I., Menzel, A., et al.
Shifting plant phenology in response to global change [J]. Trends in Ecology & Evolution,
2007, 22(7): 357–365.
[5]
Nagai, S.,
Ichie, T., Yoneyama, A., et al.
Usability of time-lapse digital camera images to detect characteristics of tree
phenology in a tropical rainforest [J]. Ecological
Informatics, 2016, 32: 91–106.
[6]
Vrieling, A.,
Meroni, M., Darvishzadeh, R., et al.
Vegetation phenology from Sentinel-2 and field cameras for a Dutch barrier
island [J]. Remote Sensing of Environment,
2018, 215(15): 517–529.
[7]
Zhang, X.,
Jayavelu, S., Liu, L., et al.
Evaluation of land surface phenology from VIIRS data using time series of
PhenoCam imagery [J]. Agricultural and
Forest Meteorology, 2018, 256–257: 137–149.
[8]
马勇刚,
刘素红. 中亚地区陆表物候逐年数据集(1982–2015) [DB/OL]. 全球变化数据仓储, 2020. DOI: 10.3974/geodb.2020.01.05.V1.
[9]
全球变化科学研究数据出版系统. 全球变化科学研究数据共享政策[OL]. DOI: 10.3974/dp.policy. 2014.05 (2017年更新).
[10]
Jonsson, P.,
Eklundh, L. Seasonality extraction by function fitting to time-series of
satellite sensor data [J]. IEEE
Transactions on Geoscience and Remote Sensing, 2002, 40(8): 1824–1832.
[11]
Zhang, X.,
Friedl, M. A., Schaaf, C. B., et al.
Monitoring vegetation phenology using MODIS [J]. Remote Sensing of Environment, 2003, 84(3): 471–475.
[12]
McLeod, A. I.
Kendall rank correlation and Mann-Kendall trend test [Z]. R Package “Kendall”, 2005.
http://www.stats.uwo.ca/faculty/aim.